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Theory (TST) limit is identical to rigorous Quantum Transition-State Theory, and we find that its long¬ 
time limit is independent of the location of the dividing surface. TRPMD rate theory is then applied to 
one-dimensional model systems, the atom-diatom bimolecular reactions H+H 2 , D+MuH and F+H 2 , and the 
prototypical polyatomic reaction H+CH 4 . Above the crossover temperature, the TRPMD rate is virtually 
invariant to the strength of the friction applied to the internal ring-polymer normal modes, and beneath the 
crossover temperature the TRPMD rate generally decreases with increasing friction, in agreement with the 
predictions of Kramers theory. We therefore find that TRPMD is approximately equal to, or less accurate 
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I. INTRODUCTION 


The accurate computation of thermal quantum rates 
is a major challenge in theoretical chemistry, as a purely 
classical description of the kinetics fails to capture zero- 
point energy, tunnelling, and phase effect^®. Exact 
solutions using correlation functions, developed by Ya¬ 
mamoto, Miller, and other^Hll 

are only tractable for 
small or model systems, as the difficulty of computation 
scales exponentially with the size of the system. 

Consequently, numerous approximate treatments have 
been developed, which can be broadly classed as those 
seeking an accurate description of the quantum statistics 
without direct calculation of the dynamics, and those 
which also seek to use an approximate quantum dy¬ 
namics. Methods in the first category include instan- 
ton theorjPEl, “quantum instanton d 18 ! 19 ! and various 
transition-state theory (TST) approached? 025 . Of many 
approximate quantum dynamics methods, particularly 
successful ones include the linearized semiclassical initial- 
value representation (LSC-IVRj®®, centroid molecular 
dynamics (CMDjPp®, and ring polymer molecular dy¬ 
namics (RPMDp^^Sl. 

RPMD has been very successful for the computation 
of thermal quantum rates in condensed-phase processes, 
due to the possibility of implementation in complex sys¬ 
tems such as (proton-coupled) electron transfer reac¬ 
tion dynamics or enzym e catalys is,®® and especially in 
small gas-phase systemd®®® where comparison with 
exact quantum rates and experimental data has demon¬ 


strated that RPMD rate theory is a consistent and re¬ 
liable approach with a high level of accuracy. These 
numerical results have shown that RPMD rate theory 
is exact in the high-temperature limit (which can also 
be shown algebraically®), reliable at intermediate tem¬ 
peratures, and more accurate than other approximate 
methods in the deep tunnelling regime (see Eq. (24) 
below), where it is within a factor of 2-3 of the exact 
quantum result. RPMD also captures zero-point energy 
effects,® and provid es very accurate estimates for barri¬ 
erless reaction d 48 ! 59 ! It has been found to systematically 
overestimate thermal rates for asymmetric reactions and 
underestimate them for symmetric (and quasisymmet- 
ric) reactions in the deep tunnelling regime (Note that 
zero-point energy effects along the reaction coordinate 
must be taken into account when assigning the reaction 
symmetry. | 13 l 52 l Recently a general code for RPMD cal¬ 
culations (RPMDrate) has been developed.® 


Another appealing feature of RPMD rate theory is 
its rigorous independence to the location of the divid¬ 
ing surface between products and reactants 8 -^, a prop¬ 
erty shared by classical rate theory and the exact quan¬ 
tum ratd®, but not by many transition-state theory ap¬ 
proaches. The t — t 0+, TST limit of RPMD (RPMD- 
TST) is identical to true QTST: the instantaneous ther¬ 
mal quantum flux through a position-space dividing sur¬ 
face which is equal to the exact quantum rate in the 
absence of recrossinj^HSSJ A corollary of this is that 
RPMD will be exact for a parabolic barrier (where there 
is no recrossing of the optimal dividing surface by RPMD 
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dynamics or quantum dynamics, and QTST is therefore 
also exact (PR 

When the centroid is used as the dividing surface 
(see Eq. (25) below), RP MD-TST reduces to the ear¬ 
lier theory of centroid-TST 25 3C ® 6 ^, which is a good 
approximation for symmetric barriers but significantly 
overestimates the rate for asymmetric barriers at low 
temperatureJ^HHHill. This effect is attributable to the 
centroid being a poor dividing surface beneath the 
‘crossover’ temperature into deep tunnelling 13 } In this 
‘deep tunnelling’ regime, RPMD-TST has a close rela¬ 
tionship to semiclassical “Im F” instanton theory^ 3 -S3 
which has been very successful for calculating rates be¬ 
neath the crossover temperature, though has no first- 
principles derivation 14 ! and was recently shown to be less 
accurate than QTST when applied to realistic multidi¬ 
mensional reaction! 70 . 


Very recently, both CMD and RPMD have been 
obtained from the exact quantum Kubo-transformecf 71 ! 
time-correlation function (with explicit error t erms ) via a 
Boltzmann-conserving “Matsubara dynamics’^ 73 which 
considers evolution of the low-frequency, smooth “Mat- 
subara” modes of the path integral 7 ^. Matsubara dynam¬ 
ics suffers from the sign problem and is not presently 
amenable to computation on large systems. However, 
by taking a mean-field approximation to the centroid 
dynamics, such that fluctuations around the centroid 
are discarded, one obtains CMD .' 3 Alternatively, if the 
momentum contour is moved into the complex plane 
in order to make the quantum Boltzmann distribution 
real, a complex Liouvillian arises, the imaginary part 
of which only affects the higher, non-centroid, normal 
modes. Discarding the imaginary Liouvillian leads to 
spurious springs in the dynamics and gives RPMD .' 3 
Consequently, RPMD will be a reasonable approximation 
to Matsubara dynamics, provided that the timescale over 
which the resultant dynamics is required (the timescale 
of ‘falling off’ the barrier in rate theory) is shorter than 
the timescale over which the springs ‘contaminate’ the 
dynamics of interest (in rate theory, this is usually cou¬ 
pling of the springs in the higher normal modes to the 
motion of the centroid dividing surface via anharmonicity 
in the potential). 

Both RPMD and CMD are inaccurate for the compu¬ 
tation of multidimensional spectra: the neglect of fluctu¬ 
ations in CMD leads to the “curvature problem” where 
the spectrum is red-shifted and broadened, whereas in 
RPMD the springs couple to the external potential lead¬ 
ing to “spurious resonances”* 75 ^ 76 : Recently, this problem 
has been solved by attaching a Langevin thermostat ' 7 
to the internal modes of the ring polymer^ (which had 
previously been used for the computation of statistical 
propertied 79 !), and the resulting Thermostatted RPMD 
(TRPMD) had neither the curvature nor resonance prob¬ 
lem. 


The success of RPMD for rate calculation, and the 
attachment of a thermostat for improving its compu¬ 
tation of spectra, naturally motivates studying whether 


TRPMD will be superior for the computation of ther¬ 
mal qu antum rates to RPMD (and other approximate 
theories)P 117 ’^, which this article investigates. Given that 
RPMD is one of the most accurate approximate meth¬ 
ods for systems where the quantum rates are available 
for comparison, further improvements would be of con¬ 
siderable benefit to the field. 


We firstly review TRPMD dynamics in section IIA 


followed by developing TRPMD rate theory in sec¬ 
tion |II B| To predict the behaviour of the RPMD rate 
compared to the TRPMD rate, we apply one-dimensional 
Kramers theory 861 to the ring-polymer potential energy 
surface in section III Cl Numerical results in section m 
apply TRPMD to the symmetric and asymmetric Eckart 
barriers followed by representative bimolecular reac¬ 
tions: H+H 2 (symmetric), D+MuH (quasisymmetrical), 
H+CH 4 (prototypical polyatomic reaction) and F+H 2 
(asymmetric and highly anharmonic). Conclusions and 


avenues for further research are presented in section IV 


II. THEORY 

A. Thermostatted Ring Polymer Molecular Dynamics 

For simplicity we consider a one-dimensional system 
(F = 1 ) with position q and associated momentum p 
at inverse temperature (3 = l/feeT, where the iV-bead 
ring-polymer Hamiltonian i j SS I Sl I I 

N—l 2 

H N ( p,q) = E^ + ^(q) (!) 

0 

with the ring-polymer potential 


N-l 

Un( q) = ^2 \mu 2 N (q t - q,- 1) 2 + V(qi) (2) 

i =0 

and the frequency of the ring-polymer springs lon = 
1/PnH, where (3n = /3/N. Generalization to further di¬ 
mensions follows immediately, and merely requires more 
indices.^ 

The ring polymer is time-evolved by pr opaga ting 
stochastic trajectories using TRPMD dynamic d 78 * 79 } 


P = - VqtMq) 



r P + 



( 3 ) 

( 4 ) 


where q = (qo,..., qN-i) is the vector of bead positions 
and p the vector of bead momenta, with V q the grad 
operator in position-space, £(f) a vector of N uniform 
Gaussian deviates with zero mean and unit variance, and 
r the N x N positive semi-definite friction matrix 78 } 
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The Fokker-Planck operator corresponding to the 
TRPMD dynamics in Eqs. Q and Q ip2 

An = — — • V q + f/iv(q)V q ■ \^ p 

TO 

TTL 

+ Vp . r . p+ _ v p . r .V p (5) 

(where the arrows correspond to the direction in which 
the derivative act^~) and for any T, TRPMD dynam¬ 
ics will conserve the quantum Boltzmann distribution 
(A N e~^ HN ^ = 0), a feature shared by RPMD and 
CMD but not som e other approximate methods such as 
LSC- 1 VFpmZEUIU. We then show in appendix [A] that 
TRPMD obeys detailed balance, such that the TRPMD 
correlation function is invariant to swapping the opera¬ 
tors at zero time and finite time, and changing the sign 
of the momenta. 

The time-evolu tion of an observable is given by the 
adjoint of Eq. 

A^n = q Lbv(q)^q ' 

TO 

vn 

~ p ' r ' Vp + ^ Vp ' r ' Vp ' (6) 

In the zero-friction limit, T = 0 and _4.jy = C' N , where 
£^ is the adjoint of the Liouvillian corresponding to 
deterministic ring-polymer trajectories 73 . 


B. TRPMD rate theory 

We assume the standard depiction of rate dynamics, 
with a thermal distribution of reactants and a dividing 
surface in position space. In what follows we assume 
scattering dynamics, with the potential tending to a con¬ 
stant value at large separation of products and reactants. 
The methodology is then immediately applicable to con¬ 
densed phase systems subject to the usual caveat that 
there is sufficient sepa ratio n of timescales between reac¬ 
tion and equilibrationP^M) 

The exact quantum rate can be formally given as 
the long-time limit of the flux-side time-correlation 
functioffS^ 


£qm(/3) = 


„qm 


lim 

£->-oo 


it) 


Q*iP) 


(7) 


where QA§) is the partition function in the reactant re¬ 
gion andS3J 


„QM 

C fs 


it) 


1 

P 



da Tr 


e ~(P~a)H p e ~aH ^Ht/h^-iHt/h 


( 8 ) 


with F and h. the quantum flux and side operators re¬ 
spectively, and H the Hamiltonian for the system. The 
quantum rate can equivalently be given as minus the 


long-time limit of the time-derivative of the side-side cor¬ 
relation function, or the integral over the flux-flux corre¬ 
lation functiorP. 

The TRPMD side-side correlation function is 

c ™ D( ‘>=(2 •k*J d °h 

x e~^ H ^h[f(q)]h[f(q t )} (9) 

where / dq = dq 0 dqi ... dq N - l and like¬ 
wise for f dp, and q t = q t (p, q, t) is obtained by evolu¬ 
tion of (p, q) for time t with TRPMD dynamics. The ring 
polymer reaction co-ordinate /(q) is defined such that 
the dividing surface is at /(q) = 0, and that /(q) > 0 
corresponds to products and /(q) < 0 to reactants. 

Direct differentiation of the side-side correlation func¬ 
tion using the Fokker-Planck operator in Eq. © yields 
the TRPMD flux-side time-correlation function 


CfeRPMD^) = 

- ^cl nPMD (t) 
dt ss w 

( 10 ) 

1 

J dp j dq e -A^(p,q) 


(2 ttH) n 



x d[f{q)}S N (p,q)h[f(q t )) 

( 11 ) 


where Sn( P, q) is the flux perpendicular to /(q) at time 
t = 0, 


JV-l 


Sn( p,q) = 


2 — 0 


dfj q) Vi 
dqt m 


( 12 ) 


We approximate the long-time limit of the quantum flux- 
side time-correlation function in Eq. ([ 8 ]) as the long-time 
limit of the TRPMD flux-side time-correlation function 


in Eq. (11), leading to the TRPMD approximation to the 


quantum rate as 


/^TRPMD/i’i 

&TRPMD iP) = lim fs , ' -• 

t->oo Qi(p) 


(13) 


The flux-side time-correlation function Eq. (Ill will 


decay from an initial TST (t —> 0+) value to a plateau, 
which (for a gas-phase scattering reaction with no friction 
on motion out of the reactant or product channel) will 
extend to infinity. For condensed-phase reactions (and 
gas-phase reactions with friction in exit channels) a rate 
is defined provided that there is sufficient separation of 
timescales between reaction and equilibration to define a 
plateau in (TTRPMD ^ [83] w ] 1 j c ] :l a t ver y ] on g times (of the 
order fcxRPMD(/3 )~ 1 for a unimolecular reaction) tends to 
zerc jll] 

Further differentiation of the flux-side time-correlation 
function (with the adjoint of the Fokker-Planck opera¬ 
tor in Eq. Q) yields the TRPMD flux-flux correlation 
function 


/"rTRPMD 

Off 


(*) = 


1 


dp dq e 


-0n ffjv(pq) 


( 2 nh) N j j 

x £[/(q)]SW(p, <i)6[f(<h)]SN (p*,qt) 

(14) 
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which, by construction, must be zero in the plateau re¬ 
gion, during which no trajectories recross the dividing 
surface. 

Like RPMD rate theory, TRPMD has the appealing 
feature that its short-time (TST) limit is identical to true 
Quantum Transition-State Theory (QTST), as can be 
observed by applying the short-time limit of the Fokker- 
Planck propagator e- 4 "* to /(q), yielding® 


lim 

t —>-0-)- 


C TRPMD( t ) 


Qt(P) 



(15) 


2Sij sm(jn/N)/ftiyh, the friction matrix is given by 

r = 2ATJ1T t . (17) 

Here A is an adjustable parameter, with A = 1 giving crit¬ 
ical damping of the free ring polymer vibrations, A = 0.5 
corresponding to optimal sampling of the free ring poly¬ 
mer potential energy, and A = 0 corresponding to zero 
friction (i.e. RPMD).' 8 ,n A crucial consequence of this 
choice of friction matrix is that the centroid of the ring 
polymer is unthermostatted, and the short-time error 
of TRPMD from exact q uantu m dynamics is therefore 
0(t 7 ), the same as RPMD.' 8 


where &q M (/3) is the QTST rat^2JMEH_ i n Appendix |b 


V QM 

we then show that the TRPMD rate in Eq. (Il3j) is rig¬ 
orously independent of the location of the dividing sur¬ 
face. Consequently, the TRPMD rate will equal the exact 
quantum rate in the absence of recrossing of the opti¬ 
mal dividing surface (and those orthogonal to it in path- 
integral space) by either the exact quantum or TRPMD 
dynamics.^ We also note that Eq. (151 holds regardless 
of the value of the friction matrix T and that recrossing 
of individual (stochastic) trajectories can only reduce the 
TRPMD rate from the QTST value, and hence QTST is 
an upper bound to the long-time TRPMD rate. 


In the following calculations we use a friction matrix 
which corresponds to damping of the free ring polymer 
vibrational frequencies, and which has b een used in previ¬ 
ous studies of TRPMD for spectraP^HIl For an orthogonal 
transformation matrix T such that 


T t KT = mn 2 (16) 


where K is the spring matrix in Eq. ([2]) and fly = 


C. Relation to Kramers Theory 

To provide a qualitative description of the effect of 
friction on the TRPMD transmission coefficient, we ap¬ 
ply classical Kramers theorjEO in the extended NF- 
dimensional ring polymer space, governed by dynamics 
on the (temperature-dependent) ring-polymer potential 
energy surface in Eq. <§• Since the short-time limit of 
TRPMD rate theory is equal to QTST, and its long-time 
limit invariant to the location of the dividing surface, 
TRPMD will give the QTST rate through the optimal 
dividing surface (defined as the surface which minimises 
^qm(i®)W weighted by any recrossings of that dividing 
surface by the respective dynamics. We express this us¬ 
ing the Bennett-Chandler factorizatiorPil, 

fcTRPMD(/3) = ^QM(^) ft TRPMD(^) (18) 

where &q M (/ 3) is the QTST rate, the asterisk denotes 
that the optimal dividing surface /*(q) is used and the 
TRPMD transmission coefficient is 


* / rip / rfq e-P NHN ^6\r(g)]S* N (p, q)ft[/*(q t )] 

^trpmd! J fdpfdq e-^MP^S[f*( q )]S* N (p, q )h[S^(p, q )] 


(19) 


with analogous expressions to Eqs. (18) and (19) for 
RPMD. To examine the explicit effect of friction on the 
TRPMD rate we define the ratio 


saddle point on the ring-polymer potential energy sur¬ 
face, in which case the TRPMD transmission c oeffici ent 
can be approximated by the Kramers expressio d — I— :91 


/ o\ &TRPMD (/?) 
xm = « 


( 20 ) 


lim KtrpmdW 

t->-oo 



RP a RP 


( 22 ) 


and from Eq. (18) 


xm = lim K msM, 

t->°o K RPMD (f) 


( 21 ) 


We then assume that the recrossing dynamics is dom¬ 
inated by one-dimensional motion through a parabolic 


where formally a R p = 7 rp/ 2 u; R p, with y R p the fric¬ 
tion along the reaction co-ordinate and w R p the bar¬ 
rier frequency in ring-polymer space. For a general F- 
dinrensional system finding /*(q) and thereby comput¬ 
ing 7 rp and w R p is largely intractable. However, we 
expect y R p oc A, and therefore define d RP = curp/A 
where the dimensionless parameter b R p is expected to 
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be independent of A for a given system and tempera¬ 
ture, and represents the sensitivity of the TRPMD rate 
to friction. We further approximate that there is min¬ 
imal recrossing of the optimal dividing surface by the 
(unthermostatted) ring polymer trajectories such that 
lim^oo 4 pmd (t) ~ 1, 921 leading to 


X\(H) - \Jl + A 2 5|p - Aa RP . 


(23) 


Equation (23) relates the ratio of the TRPMD and 


RPMD rates as a function of A with one parameter Arp, 
and without requiring knowledge of the precise location 
of the optimal dividing surface /*(q). However, we can 
use general observations concerning which ring-polymer 
normal modes contribute to /* (q) to determine the likely 
sensitivity of the TRPMD rate to friction. Above the 
crossover temperature into deep tunnelling, defined by 13 


He = 


27r 
hu>b 


(24) 


where Ub is the barrier frequency in the external potential 
V(q), the optimal dividing surface is well approximated 
by the centroid^ 


f (q) 


1 

N 


N -1 


* 




(25) 


where is the maximum in V{q). As the centroid is not 
thermostatted (since Qqo = 0 ), in this regime 7 rp = 0 = 
Arp and we therefore predict from Eq. (23) that the rate 
will be independent of A, i.e. fcT R PMD(/oTcs ^rpmd(/?)■ 
Beneath the crossover temperature, the saddle point on 
the ring-polymer potential energy surface bends int o the 
space of the first degenerate pair of normal modesP 3 ^ 
For symmetric systems, the optimal dividing surface is 
still the centroid expression in Eq. (25) and (insofar as the 
reaction dynamics can be considered one-dimensional) 
Arp ~ 0, so & T rpmd(/3) — ^rpmd(/3)- 
For asymmetric reactions, the optimal dividing surface 
is now a function of both the centroid and first degen¬ 
erate pair of normal modes (which are thermostatted) ^ 
and we expect Arp > 0. From Eq. (23) the TRPMD rate 
will decrease linearly with A for small A, for large friction 
as A -1 , and the ratio of the TRPMD to RPMD rates to 
be a convex function of A. This behaviour would also 
be expected for symmetric reactions beneath the second 
crossover temperature where the optimal dividing sur¬ 
face bends into the space of the second degenerate pair 
of normal modest In all cases one would expect that in¬ 
creasing friction would either have no effect on the rate, 
or at sufficiently low temperatures cause it to decrease. 

It should be stressed that Eq. (23) is a considerable 
simplification of the TRPMD dynamics and is not ex¬ 
pected to be reliable in systems where the ring poly¬ 
mer potential energy surface is highly anharmonic or 
skewed (such as F+H 2 investigated below). In fact, 
even for a one-dimensional system, the minimum energy 


path on the ./V-dimensional ring polymer potential energy 
surface shows a significant skew beneath the crossover 
temperatur^H. The utility of Eq. (23) lies in its sim¬ 
plicity and qualitative description of friction-induced re¬ 
crossing. 


III. RESULTS 

We initially study the benchmark one-dimensional 
symmetric and asymmetric Eckart barriers before pro¬ 
gressing to the multidimensional reactions H+H 2 (sym¬ 
metric), D+MuH (quasisymmetrical), H+CH 4 (asym¬ 
metric, polyatomic) and F+H 2 (asymmetric, anhar¬ 
monic) . 


A. One-dimensional results 


The methodology for computation of TRPMD reac¬ 
tion rates is identical to that of RPME®, except for the 
thermostat attached to the internal normal modes of the 
ring polymer, achieved using the algorithm in Ref. 1791 
The Bennett-Chandlei® factorization was employed, and 
the same dynamics can be used for thermodymamic in¬ 
tegration along the reaction co-ordinate (to calculate the 
QTST rate) as to propa gate t rajectories (to calculate the 

transmission coefficient)! 78 * 79 * _ 

We firstly examine the symmetric Eckart barriei^lMJ, 

V(q) = V 0 sech 2 {q/a), (26) 


and to facilitate comparison with the literature 13 !®!^, 
use parameters to model the H+H 2 reaction: Vo = 
0.425eV, a = 0.734a 0 , and m = 1061m e , leading to 
a crossover temperature of ksHc = 2.69 x 10 _ 3 K _1 . 
The centroid reaction co-ordinate of Eq. (25) was used 
throughout. Results for a variety of temperatures and 
values of friction parameter A are presented in Fig. |T| 
and values of Arp obtained by nonlinear least squares in 
Table ID 

Slightly beneath the crossover temperature (fcB/3 c = 
3 x 10 _ 3 K _1 ), the TRPMD rate is indepedent of the 
value of friction (Arp = 0), as predicted by Kramers 
theory. Some sensitivity to A is seen before twice the 
crossover temperature, which is likely to be a breakdown 
of the one-dimensional assumption of Kramers theory; 
while the centroid is the optimal dividing surface, the 
minimum energy path bends into the space of the (ther¬ 
mostatted) lowest pair of normal modes***. Beneath twice 
the crossover temperature the friction parameter has a 
significant effect on the rate, as to be expected from the 
second degenerate pair of normal modes becoming part 
of the optimal dividing surfac^^D. The functional form 
°f Xx{H) is also in accordance with the predictions of 
Kramers theory, monotonically decreasing as A rises, and 
being a convex function of A. 

Since RPMD underestimates the rate for this sym¬ 
metric reaction (and many other^®), adding friction to 
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TABLE I. Dimensionless friction sensitivity parameter Srp from Eq. (231, fitted by nonlinear least squares to simulation data. 



ID Eckart barriers 


Multidimensional reactions 


fcBd/HT'TC- 1 3 

5 

7 

T/K 

500 

300 

200 

Symmetric 

<0.01 

0.11 

0.37 

h+h 2 

0.01 

0.16 

0.45 





D+MuH 

0.20 

0.45 

0.71 


d/a.u. 4 

8 

12 

T/K 

500 

300 

200 

Asymmetric 

0.00 

0.06 

0.17 

H+CH 4 

0.00 

0.10 

0.16 (250K) 





f+h 2 

-0.01 

-0.01 

0.00 


4.0E-3-, 


3.0E-3 -f. 
2.0E-3 - 
1.0E-3- 
0.0 


-x -* - 


P = 3 


•X-- X --X-X 

X TRPMD(A) 

-QM 

- Kramers 


0.00 0.25 0.50 0.75 1.00 1.25 1.50 



3.0E-6 n 
2.0E-6 


1.0E-6> x - * - * - x __ 

P = 7 * 

0.0 


0.00 0.25 0.50 0.75 1.00 1.25 1.50 


2.0 

1.5 

1.0 

0.5 

0.0 


■ 

Sr - A - 

- P = 4 

1 1 

- -x- 

1 

- -K- 

1 

-;*- X- - X 

x TRPMD(A) 

- QM 

- Kramers 

1 1 1 1 1 1 


0.00 0.25 0.50 0.75 1.00 1.25 1.50 






FIG. 1. Results for the symmetric Eckart barrier, showing the 
TRPMD result as a function of A (red crosses), fitted Kramers 
curve (green dashes) and quantum result (black line). /? is 
quoted in units of fcg 1 10 _3 K _1 and the crossover temperature 
is fe B /3 c = 2.69 x l(r 3 K _1 . 


FIG. 2. Results for the asymmetric Eckart barrier quoted 
as c(/3 ) [Eq. (281], showing the TRPMD result as a function 
of A (red crosses), fitted Kramers curve (green dashes) and 
quantum result (black line). The crossover temperature is 
f) c = 27ra.u. 


R.PMD decreases its accuracy in approximating the quan¬ 
tum rate for this system. 

The asymmetric Eckart barrier is given by 38 


V(q) 


A B 

1 + e- 2 9/“ + cosh 2 (g/a) 


(27) 


where A = —18/n, B = 13.5/7T and a = 8/\/37r in atomic 
units (h = fcs = m = 1), giving a crossover tempera¬ 
ture of = 2 tt. To facilitate comparison with previous 
literatur e! 1 ’d’ 38 l 94 l H -^ the results are presented in Fig. [ 2 ] as 


the ratio 


c(/3) 


m 

^-clas {P) 


(28) 


and cIrp values in Table |Tj 

Above the crossover temperature, TRPMD is invari¬ 
ant to the value of the friction parameter, and beneath 
the crossover temperature, increasing A results in a de¬ 
crease in the rate, such that TRPMD is closer to the 
exact quantum result than RPMD for all A > 0 in this 
system. The decrease in the TRPMD rate with A is qual¬ 
itatively described by the crude Kramers approximation 
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q* 


FIG. 3. TRPMD (green dashes), RPMD (blue dots) and 
QTST (centroid dividing surface, red line) rates for the asym¬ 
metric Eckart barrier at /3 = 12, as a function of the dividing 
surface qA 


(see Fig. [2]), and it therefore seems that the improved ac¬ 
curacy of TRPMD could be a fortuitous cancellation be¬ 
tween the overestimation of the quantum rate by QTST, 
and the friction-induced recrossing of the optimal divid¬ 
ing surface by TRPMD trajectories. There is no particu¬ 
lar a priori reason to suppose that one value of A should 
provide superior results; from Fig. [2] at /3 = 8 a fric¬ 
tion parameter of A = 1.25 causes TRPMD to equal the 
quantum result to within graphical accuracy, whereas at 
/3 = 12 this value of friction parameter causes overesti¬ 
mation of the rate, and further calculations (not shown) 
show that A = 5 is needed for TRPMD and the quantum 
rates to agree. 


The numerical results also show a slightly higher cur¬ 
vature in A;trpmd(/3) as a function of A than Eq. (23) 
would predict, suggesting that the TRPMD rate reaches 
an asymptote at a finite value, rather than at zero as 
the Kramers model would suggest. We suspect this is 
a breakdown of one-dimensional Kramers theory, since 
in the A —> 00 limit the system can still react via the 
unthermostatted centroid co-ordinate, but may have to 
surmount a higher barrier on the ring polymer potential 
energy surface. 


We then investigate the effect of changing the loca¬ 
tion of the centroid dividing surface on the TRPMD rate. 
RPMD is already known to be invariant to the location 
of the dividing surfaces!, and we therefore choose a sys¬ 
tem for which TRPMD and RPMD are likely to differ 
the most, namely a low-temperature, asymmetric sys¬ 
tem where there is expected to be significant involvement 
of the thermostatted lowest degenerate pair of normal 
modes in crossing the barrier. The asymmetric Eckart 
barrier at /3 = 12 is therefore used as a particularly 
harsh test, with the result plotted in Fig. [3] Although 
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FIG. 4. Results for the H+FF reaction as a function of A. 
Kramers is the fitted Kramers curve (see text). The crossover 
temperature is 345K. 


the centroid-density QTST result varies by almost a fac¬ 
tor of six across the range of dividing surfaces considered 
(—3 < < —2a.u.), both the TRPMD and RPMD rates 

are invariant to the location of the dividing surface. We 
also observe that, even with the optimal dividing sur¬ 
face, centroid -dens ity QTST significantly overestimates 
the exact rat c! 38 l 95 l 


B. Multidimensional results 

The results are calculated using adapted RPMDrate 
cod^H, with details summarized in Table |n| In the cal¬ 
culations reported below we used the potential energy 
surface developed by Boothroyd et al. (BKMP2 PES) 
for H+H 2 and D+MuHpl the Stark-Werner (SW) po¬ 
tential energy surface for F+FQp^ and the PES-2008 
potential energy surface developed by Corchado et al. 
for H+CH 4 PD The computation of the free energy was 
achieved using umbrella integration 100 101 with TRPMD 
and checked against sta nda rd umbrella integration with 
an Andersen thermostats^. 

H+H 2 represents the simplest atom-diatom scatter¬ 
ing reaction and has been the subject of numerous 




















TABLE II. Input parameters for the TRPMD calculations on the H + H 2 , D + MuH, and F + H 2 reactions. The explanation of the format of 
the input file can be found in the RPMDrate code manual (see Ref. 82 and http://www.mit.edu/ ysuleyma/rpmdrate). 


Parameter 



Reaction 


Explanation 


H + H 2 

D + MuH 

F + H 2 

H + CH 4 


Command line parameters 

Temp 

200; 300; 500 


250; 300; 500 

Temperature (K) 

Nbeads 

128 

512 

384 (200 K) 

192 (250 K) 

Number of beads in the TRPMD calculations 




256 (300 K) 

128 (300 K) 





64 (500 K) 

64 (500 K) 


Dividing surface parameters 

Ro 0 

30 

30 

30 

30 

Dividing surface si parameter (ao) 

Abonds 

1 

1 

1 

1 

Number of forming and breaking bonds 

A”channel 

2 

1 

2 

4 

Number of equivalent product channels 

Thermostat options 

thermostat 

’GLE/Andersen’ 



Thermostat for the QTST calculations 

A 

0; 0.25; 

0.5; 0.75; 1.0; 1.5 



Friction coefficient for the recrossing factor calculations 

Biased sampling parameters 

windows 

111 

111 

111 

111 

Number of windows 


-0.05 

-0.05 

-0.05 

-0.05 

Center of the first window 

dt 

0.01 

0.01 

0.01 

0.01 

Window spacing step 

Cjv 

1.05 

1.05 

1.05 

1.05 

Center of the last window 

dt 

0.0001 

0.0001 

0.0001 

0.0001 

Time step (ps) 

ki 

2.72 

2.72 

2.72 

2.72 

Umbrella force constant ((T/K) eV) 

N trajectory 

200 

200 

200 

200 

Number of trajectories 

^equilibration 

20 

20 

20 

20 

Equilibration period (ps) 

^sampling 

100 

100 

100 

100 

Sampling period in each trajectory (ps) 

Ni 

2 x 10 8 

2 x 10 8 

2 x 10 8 

2 x 10 s 

Total number of sampling points 

Potential of 

mean force calculation 




£0 

- 0.02 

- 0.02 

- 0.02 

- 0.02 

Start of umbrella integration 

** 

1 . 0000 * 

0.9912 (200 K)* 

0.9671 (200 K)* 

1.0093 (250 K)* 

End of umbrella integration 



0.9904 (300 K)* 

0.9885 (300 K)* 

1.0074 (300 K)* 




0.9837 (500 K)* 

0.9947 (500 K)* 

1.0026 (500 K)* 


Abins 

5000 

5000 

5000 

5000 

Number of bins 

Recrossing factor calculation 

dt 

0.0001 

0.00003 

0.0001 

0.0001 

Time step (ps) 

^equilibration 

20 

20 

20 

20 

Equilibration period (ps) in the constrained (parent) 






trajectory 

Atotalchild 

100000 

100000 

500000 

500000 

Total number of unconstrained (child) trajectories 

tchildsampling 

20 

20 

20 

20 

Sampling increment along the parent trajectory (ps) 

Achild 

100 

100 

100 

100 

Number of child trajectories per one 






initially constrained configuration 

^child 

0.05 

0.2 

0.2 

0.1 

Length of child trajectories (ps) 


* Detected automatically by RPMDrate. 


studie d 44152 1 11 ^ 1 . The PES is symmetric and with a rel¬ 
atively large skew angle (60°), and a crossover tempera¬ 
ture of 345K. The results in Fig. [3] show that the rate is 
essentially invariant to the value of A above the crossover 
temperature. At 300K there is a slight decrease in the 
rate with increasing friction from 0 to 1.5 (~25 %), and 
this is far more pronounced at 200K where the A = 1.5 
result is almost half that of the A = 0 (RPMD) result. 

D+MuH is “quasisymmetrical” since DMu and MuH 
have very similar zero-point energies, and one would 
therefore expect the RPMD rate to underestimate the ex¬ 
act quantum rate^. Since it is Mu-transfer the crossover 
temperature is very high (860 K) and therefore this reac¬ 


tion can be considered as a stress test for the deep tun¬ 
neling regime. The results in Fig. [5] show that friction in 
the TRPMD dynamics causes further underestimation of 
the rate, especially at low temperatures; for A = 1.5 at 
200K, TRPMD underestimates the exact quantum rate 
by almost an order of magnitude, and even at 500K it 
decreases by ~40% over the range of A explored here. 

As an example of a typical asymmetric reaction, re¬ 
sults for H+CH 4 are plotted in Fig. [ 6 j which has a 
crossover temperature of 341K. RPMD is well-known to 
overestimate the quantum rate for this system at low 
temperatures PS! Fig. [ 6 ] shows that above the crossover 
temperature (500K) the friction parameter has a negligi- 
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FIG. 5. As for Fig. [4] but for the D+MuH reaction with a FIG. 6. Results for the H+CH4 reaction. The crossover tem- 
high crossover temperature of 860K. perature is 341K. 


ble effect on the rate. As the temperature is decreased 
below the crossover temperature (300K and 250K), the 
friction induces more recrossings of the dividing surface 
and, as a result, the TRPMD rate approaches the exact 
quantum rate with increasing the friction parameter. 

Thus far, Kramers theory has been surprisingly suc¬ 
cessful at qualitatively explaining the behaviour of the 
TRPMD rate with increasing friction. Present results 
would suggest that TRPMD would therefore improve 
upon RPMD for all asymmetric reactions, where RPMD 
generally overestimates the rate beneath crossove j 13 l 52 l 
We then examine another prototypical asymmetric reac¬ 
tion, F+H 2 , with a low crossover temperature of 264K. 
Fig. [7] shows that at 500K and 300K, the TRPMD rate 
is in good agreement with the quantum result, but in¬ 
creases very slightly with A causing a spurious small neg¬ 
ative value of 5 rp in Table [ij Beneath crossover, at 200K 
the rate is virtually independent of lambda, apart from 
a very slight increase around A = 0.5. Consequently, 
TRPMD fares no better than RPMD for this system, 
contrary to the H+CH 4 results and the predictions of 
Kramers theory. This is likely attributable to a highly 
anharmonic and exothermic energy p rofile, and a very 
flat saddle point in ring-polymer spac^SIifiil. 

As can be seen from the graphs, the simple Kramers 


prediction is in surprisingly good qualitative agreement 
with the numerical results (apart from F+H 2 beneath 
crossover), even for the multidimensional cases, which is 
probably attributable to those reactions being dominated 
by a significant thermal barrier which appears parabolic 
on the ring-polymer potential energy surface, meaning 
that the one-dimensional Kramers model is adequate for 
capturing the friction-induced recrossing. In Table [1] the 
<Srp values, fitted to the numerical data, show that for a 
given reaction 5 rp ~ 0 above the crossover temperature, 
and beneath the crossover temperature <5 rp increases as 
the temperature is decreased. This can be qualitatively 
explained as the optimal dividing surface becoming more 
dependent on the thermostatted higher normal modes 
as the temperature is lowered 1 ^ Not surprisingly, the 
highest value of 5 rp is observed for the highly quantum 
mechanical D+MuH reaction at 200K with cIrp = 0.71. 
This is beneath one quarter of the crossover temperature, 
and one would therefore expect that friction would have 
a very significant effect on the rate. 
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FIG. 7. Results for the anharmonic and asymmetric F+Ha 
reaction with a crossover temperature of 264K. 


IV. CONCLUSIONS 


In this paper we have, for the first time, applied Ther- 
mostatted Ring Polymer Molecular Dynamics (TRPMD) 
to reaction rate theory. Regardless of the applied fric¬ 
tion, the long-time limit of the TRPMD flux-side time- 
correlation function (and therefore the TRPMD rate) is 
independent of the location of the dividing surface , and 
its short-time limit is equal to rigorous QTST 63 ^H 78 l 
In section IIC we use Kramers theory^ to predict 
that, above the crossover temperature, the RPMD and 
TRPMD rates will be similar, and beneath crossover the 
TRPMD rate for asymmetric systems will decrease with 
A, and the same effect should be observed for symmetric 
systems beneath half the crossover temperature. 


TRPMD rate theory has then been applied to the 
standard one-dimensional model systems of the symmet¬ 
ric and asymmetric Eckart barriers, followed by the bi- 
molecular reactions H+H 2 , D+MuH, H+CH 4 and F+H 2 . 
For all reactions considered, above the crossover tem¬ 
perature the TRPMD rate is virtually invariant to the 
value of A and therefore almost equal to RPMD, as pre¬ 
dicted by Kramers theory. Beneath the crossover temper¬ 
ature, most asymmetric reactions show a decrease in the 


TRPMD rate as A is increased, and in qualitative agree¬ 
ment with the Kramers prediction in Eq. (23). A simi¬ 
lar trend is observed for symmetric reactions, which also 
show some diminution in the rate with increasing friction 
above half the crossover temperature (/3 C < /3 < 2 j3 c ), 
probably due to the skewed ring-polymer PES caus¬ 
ing a breakdown in the one-dimensional assumption of 
Kramers theory. For the asymmetric and anharmonic 
case of F+H 2 , beneath the crossover temperature there 
is no significant decrease in the rate with increased fric¬ 
tion, illustrating the limitations of Kramers theory. 

These results mean that beneath the crossover tem¬ 
perature TRPMD will be a worse approximation to the 
quantum result than RPMD for symmetric and qua- 
sisymmetrical systems (where RPMD underestimates the 
rat e! 13 * 52 *) , and TRPMD will be closer to the quantum 
rate for asymmetric potentials (where RPMD overesti¬ 
mates the rate). However, the apparent increase in accu¬ 
racy for asymmetric systems appears to be a cancellation 
of errors from the overestimation of the quantum rate by 
RPMD which is then decreased by the friction in the 
non-centroid normal modes of TRPMD, and there is no 
a priori reason to suppose that one effect should equal 
the other for any given value of A. 

Although the above results do not advocate the use 
of TRPMD rate theory as generally being more accu¬ 
rate than RPMD, TRPMD rate calculation above the 
crossover temperature may be computationally advanta¬ 
geous in complex systems due to more efficient sampling 
of the ring-polymer phase s pace by TRPMD trajecto¬ 
ries than RPMD trajectories!^. TRPMD may there¬ 
fore provide the same accuracy as RPMD rate calcula¬ 
tion at a lower computational cost, and testing this in 
high-dimensional systems where RPMD has been suc¬ 
cessful, such as complex-forming reaction, r ! 48 l 55 l 59 l 6 -d, sur¬ 
face dynamic^ 47 } and enzyme catalysitPS would be a use¬ 
ful avenue of future research. 

Future work could also include non-adiabatic 
system jS H 43 l 106 tfi l0( applying a thermostat to the 
centroid to model a bath systeirPTl, and generaliza tion s 
to non-Markovian friction using Grote-Hynes theory 111 . 

In closing, present results suggest that TRPMD can 
be used above the crossover temperature for thermally 
activated reactions, and beneath crossover further testing 
is required to assess its utility for asymmetric systems. 
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Appendix A: Detailed Balance 


For a homogeneous Markov process such as TRPMD 
for which n egat ive time is not defined^!, detailed balance 
is defined as 112 


^(p':q'^|p,q, 0)p s (p,q) 

= ^(-P, q< t\ - p', q', 0)p s (p', q') (Al) 


where p s (p,q) = e -/ 3 wi?w( p , q ) j s thg stationary distribu¬ 
tion and V{p', q', t|p, q, 0) is the conditional probability 
that a ring polymer will be found at point (p', q') at time 
t , given that is was at (p, q) at time t = 0. 

To demonstrate that Eq. a is statisfied, we rewrite 
the Fokker-Planck operator Eq. |5]) as 


JV-l 


A " = -Y. ( + 


3=0 

N-1N-1 

+ EE 


d d 


2 ^ ' dpj dPi 

j=0 j ’=o ^ w 


C{p,q)j 


(A2) 


where the vectors a(p, q) = p/m, b(p, q) = 


— £//v(q)V q — r p and the matrix C(p , q) = 2mT/ /3n- 
Note that the derivatives in Eq. (A2) act on a(p, q), 
b(p.q) or C(p,q) and whatever follows them which is 
acted upon by An- 

The necessary and sufficient conditions for detailed 


balance [Eq. (Al)] to hold, in addition to p, 
ing a stationary distribution, are then given b; 


q) be- 


a(-p, q)p s (p, q) = - a(p, q)p s (p, q) (A3) 

—b(—p, q) T p s (p, q) = - b(p, q) T p s (p, q) 

+ V p ■ C(p,q)p s (p,q) (A4) 

C(—P,q) =C(p,q) (A5) 


Condition Eq. ( |A3| ) is trivially satisfied. Provided that 
the friction matrix is even w.r.t. momenta (satisfied here 
as T is not a function of p) Eq. (A5) will be satisfied. 
Eq. (A4) becomes 


m 


(r-p) p s (p,q) =- —Vp ■ Tp s (p,q) 


(A 6 ) 


which is satisfied with p s (p,q) = e _/3jvffjv ^ p q ^ and the 
friction matrix used here. 

Given that Eq. (Al) is satisfied, for an arbitrary cor¬ 
relation function one can then show 


/-iTRPMD 


(*) 


( 2 nh)N J dp J dq J dp ' j dql e ^ JvJ?Jv(p,q) A(p, q)'P(p / , q', t|p, q, 0)-B(p', q') 
(2nh) N j dp J dq j dpl J dq ' e _feffjv(p ' q) A(-p'. q')P(p', q', t|p, q, 0)B(-p, q) 


(A7) 

(AS) 


and for the Langevin trajectories considered here, which Appendix B: Independence of &trpmd(/1) to the dividing 
are continuous but not differentiable, this means surface location 


riTRPMD 
°AB 


(t) 


dp / dq e 


-/ 3 ivffw(p,q) 


( 2 irh) N j j 

x A(p,q)i?(p t ,q ( ) 


(A9) 


1 


dp / dq e 




( 2 nh) N j j 

X A(-p t ,q t )B(-p,q) (A10) 


We use a similar methodology to that which Craig and 
Manolopoulos employed for RPME^S, and give the main 
steps here. We firstly differentiate the side-side correla¬ 
tion function in Eq. ([9]) w.r.t. the location of the dividing 
surface (or any other parameter specifying the nature 
of the dividing surface), giving 


where q t = q t (p,q, t) is the vector of positions stochas¬ 
tically time-evolved according to Eqs. ([3| and Q, and 


B(p t ,q t ) = J dp' j dq' V{p', q', i|p, q, 0)B(p', q') 

(All) 


with A(pi,q t ) similarly defined. 


^ C “W=(2 ^ S ^ 1 ^"""^ 

x R/(q)M/(q*)] + M/(q)M/(q *)]}. 

(Bl) 


Since TRPMD dynamics obeys detailed balance (as 
shown in appendix 0, and the dividing surface is only 
a function of position, the second term on the RHS of 
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Eq. (B1) is identical to the first 
2 


dqt Css (*) ~{2wh) N 


dp J dq e -0N H N(p,q.) 

df( q) 


dqt 


-S[f(q)]h[f(q t )}. (B2) 


Differentiation of Eq. ( |B2| ) w.r.t. time using Eq. 
and relating the side-side and flux-side functions using 
Eq. pi), yields 


jLr (t) = - ^ 

dqt ts[> (2 nh) N 


dp dq e 


-/3jvffjv(p,q) 


x d l^ S[f(q)]5[f(q t )}S N (p t , q t )- 


dqt 


(B3) 


Equation (B3) corresponds to a trajectory commencing 
at the dividing surface at time zero and returning to it 
at time t with non-zero velocity <Sjv(Pt, qt)- At finite 
times while there is recrossing of the barrier, there will 
be trajectories satisfying these conditions, but after the 
plateau time when no trajectories recross the barrier [cf. 
Eq. @], these conditions are clearly not satisfied, and 
the rate will be independent of the location of the divid¬ 
ing surface.® 

This proof is valid for any friction matrix which satis¬ 
fies the detailed balance conditions of appendix [XJ and 
does not require the presence of ring-polymer springs in 
the potential, so is valid for any classical-like reaction 
rate calculation using Langevin dynamics. 
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